#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit.h>

//环境：apt install libgsl-dev
//编译：gcc -o run_polynomial_fitting_gsl polynomial_fitting_gsl.c -lgsl -lgslcblas -lm
//编译：gcc -o run_polynomial_fitting_gsl polynomial_fitting_gsl.c -lgsl -lgslcblas -lm -Wall -std=c99


// 定义多项式函数
double polynomial(double x, const gsl_vector *c, size_t n) {
    double result = 0;
    for (size_t i = 0; i <= n; ++i) {
        result += gsl_vector_get(c, i) * pow(x, (int)i);
    }
    return result;
}

// 多项式拟合函数
int polynomial_fit(double *x, double *y, size_t n, size_t p, gsl_vector *c) {
    gsl_matrix *X = gsl_matrix_alloc(n, p + 1);
    gsl_vector *yvec = gsl_vector_alloc(n);
    gsl_vector *w = gsl_vector_alloc(n);
    gsl_matrix *cov = gsl_matrix_alloc(p + 1, p + 1);

    // 初始化权重向量
    for (size_t i = 0; i < n; ++i) {
        gsl_vector_set(w, i, 1.0); // 假设所有数据点的权重相同
    }

    // 填充矩阵 X 和向量 yvec
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j <= p; ++j) {
            gsl_matrix_set(X, i, j, pow(x[i], j));
        }
        gsl_vector_set(yvec, i, y[i]);
    }

    // 执行拟合
    gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n, p + 1);
    double chi2;
    gsl_multifit_wlinear(X, w, yvec, c, cov, &chi2, work);
    gsl_multifit_linear_free(work);

    // 清理内存
    gsl_matrix_free(X);
    gsl_vector_free(yvec);
    gsl_vector_free(w);
    gsl_matrix_free(cov);

    return 0; // 返回成功状态
}

int main() {
    // 示例数据点
    size_t n = 3; // 数据点的数量
    double x[] = {40979.0, 59455.0, 78519.0}; // x 值
    double y[] = {10.031000, 14.349000, 19.523001}; // y 值

    // 多项式的次数
    size_t p = 2; // 例如，拟合二次多项式

    // 存储系数的向量
    gsl_vector *c = gsl_vector_alloc(p + 1);

    // 执行拟合
    polynomial_fit(x, y, n, p, c);

    // 输出结果
    printf("拟合系数: \n");
    for (size_t i = 0; i <= p; ++i) {
        printf("c%ld = %g\n", i, gsl_vector_get(c, i));
    }

    // 使用拟合多项式计算一些值
    printf("\n拟合多项式: \n");
    for(double xi = 0.0; xi <= 3.5; xi += 0.5) {
        double yi = polynomial(xi, c, p);
        printf("f(%.1f) = %.5f\n", xi, yi);
    }

    // 清理内存
    gsl_vector_free(c);

    return 0;
}
